// STL
#include <cmath>                                  // std::round()
#include <functional>                             // std::function<>
#include <vector>                                 // std::vector<>

// jScience
#include "jrandnum.hpp"                           // rand_num_uniform_Mersenne_twister()

// this
#include "statsxx/optimization/EA/EAParam.hpp"    // EAParam
#include "statsxx/optimization/EA/evolution.hpp"  // get_elite_individuals(), mutate(), crossover_kpoint()
#include "statsxx/optimization/EA/Individual.hpp" // Individual
#include "statsxx/optimization/EA/Population.hpp" // Population
#include "statsxx/optimization/EA/utility.hpp"    // tournament_select()


//
// DESC: Generate a new population via crossover and mutation.
//
inline Population reproduce(
                            const Population                                                    &population,
                            // -----
                            const std::vector<double>                                           &f,
                            const std::vector<double>                                           &v,
                            // =====
                            const std::function<double(const int, const std::vector<double> &)> &new_gene,
                            const std::function<double(const int, const std::vector<double> &)> &mutate_gene,
                            // =====
                            const EAParam                                                       &param
                            )
{
    Population offspring;

    //=========================================================
    // INITIALIZATION
    //=========================================================

    int noffspring = population.individuals.size();

    //=========================================================
    // ELITE INDIVIDUALS
    //=========================================================

    if( param.frac_elite > 0. )
    {
        bool feasible;

        // GET FEASIBLE INDIVIDUALS
        int nelite_feasible = static_cast<int>(std::round(param.frac_elite_feasible*param.frac_elite*static_cast<double>(param.pop_size)));
        // -----
        feasible = true;

        std::vector<Individual> elite_feasible = get_elite_individuals(
                                                                       nelite_feasible,
                                                                       // -----
                                                                       feasible,
                                                                       // =====
                                                                       population,
                                                                       // -----
                                                                       f,
                                                                       v
                                                                       );

        // GET INFEASIBLE INDIVIDUALS
        //
        // NOTE: If (not enough) feasible solutions are (yet) found, supplement (for now) with infeasible ones.
        //
        // NOTE: This is consistent with B. Tessema (thesis) (end of Section 4.5).
        //
        int nelite_infeasible = static_cast<int>(std::round(param.frac_elite*static_cast<double>(param.pop_size))) - elite_feasible.size();
        // -----
        feasible = false;

        std::vector<Individual> elite_infeasible = get_elite_individuals(
                                                                         nelite_infeasible,
                                                                         // -----
                                                                         feasible,
                                                                         // =====
                                                                         population,
                                                                         // -----
                                                                         f,
                                                                         v
                                                                         );

        // ADD TO OFFSPRING
        for( auto &e : elite_feasible )
        {
            offspring.individuals.push_back(e);
        }

        for( auto &e : elite_infeasible )
        {
            offspring.individuals.push_back(e);
        }

        noffspring -= (elite_feasible.size() + elite_infeasible.size());
    }

    //=========================================================
    // CROSSOVER AND/OR MUTATION
    //=========================================================

    while( noffspring > 0 )
    {
        double p = rand_num_uniform_Mersenne_twister(0., 1.);

        if( (p < param.prob_mutate_only) || (noffspring == 1) )
        {
            Individual child = tournament_select(population, offspring, param);

            mutate(
                   child,
                   // -----
                   new_gene,
                   mutate_gene,
                   // -----
                   param
                   );

            offspring.individuals.push_back(child);

            --noffspring;
        }
        else
        {
            Individual child1, child2;

            Individual mom = tournament_select(population, offspring, param);
            Individual dad = tournament_select(population, offspring, param);

            crossover_kpoint(mom, dad, child1, child2, param);

            double pm1 = rand_num_uniform_Mersenne_twister(0., 1.);
            double pm2 = rand_num_uniform_Mersenne_twister(0., 1.);

            if( pm1 < param.prob_mutate )
            {
                mutate(
                       child1,
                       // -----
                       new_gene,
                       mutate_gene,
                       // -----
                       param
                       );
            }

            if( pm2 < param.prob_mutate )
            {
                mutate(
                       child2,
                       // -----
                       new_gene,
                       mutate_gene,
                       // -----
                       param
                       );
            }

            offspring.individuals.push_back(child1);
            offspring.individuals.push_back(child2);

            noffspring -= 2;
        }
    }

    return offspring;
}
